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SUMMARY 


Surface motion of compliant walls in drag reduction experiments has been 
analyzed. Critical comparison is made between the dynamic motion of the structure 
and the postulated mechanism of drag reduction (ref. 1) . The spectrum of surface 
motion indicates that membranes over deep cavities respond at low frequencies 
and large wavelengths. The membrane over a deep cavity is therefore found not 
to yield the desired response predicted by the postulated mechanism. The mem- 
brane over a thin air gap is found to act as a wavelength chopper, and analysis 
of the nonlinear response of that compliant surface indicates its possible 
suitability for compliant wall experiments. Periodic structures are found to 
lock in the desired wavelengths of motion, and it is found that at least in 
Kramer's (ref. 2) initial experiments they could have produced high frequency 
surface motions. Laminated structures are found to be very ineffective as com- 
pliant models, except when there is no bonding between the membrane and the 
backing. Computer programs developed for these analyses are documented in this 
report . 


INTRODUCTION 

The impetus for the compliant wall drag reduction program at NASA-Langley 
Research Center came mainly from the fact that reduction in turbulent, skin friction 
drag reported in the literature (refs. 2 to 5) may be translated into potentially 
large savings in energy for CTOL aircraft. The group effort at Langley was 
initiated with the goals of: 

1. Understanding the conditions, if any, for favorable compliant motions, 

2. Duplicating previous data under rigorously controlled test conditions, and 


3. Designing improved experiments to further understanding of the compliant 
wall drag reduction phenomenon. 

The present work under NSG. 1-236 emanated from a desire to evaluate the struc- 
tural response of compliant walls in successful drag reduction experiments . 

Certain commercial materials are identified in this paper in order to specify 
adequately which materials were investigated. in the research effort. In no 
case .does such identification imply recommendation or endorsement of the product 
by NASA, nor does it imply that the materials are necessarily the only ones or 
the best ones available for the purpose. 


THE PROPOSED FLUID-STRUCTURE INTERACTION MODEL 


Bushnell (ref. 1) has proposed a compliant wall interaction model by which 
a fully turbulent flow over a vibrating compliant wall could produce reduced skin 
friction on the wall. In the past few years, researchers have developed an under- 
standing of some details of the streaky flow adjacent to a rigid wall. The pro- 
posed mechanism of Bushnell (ref. 1) suggests that the induced aerodynamic 
pressure field due to the structural motion could impose a high frequency stab- 
ilizing modulation of the "preburst" flow and thereby reducing the rate of pro- 
duction of turbulent shear. Numerical experiments by Orszag (ref. 6) and 
unpublished data by Kendall and Collins suggest that such a mechanism may well 
be valid. 

Based on this model, the favorable structural motion has a wavelength on 

the order of X + = 100, in wall units, [A = A + v/ (/C^/2 U^) ; where Cp is the 

skin friction coefficient, v is the kinematic viscosity, and U^ is the mean 

~ + 

flow velocity], a wave speed C = 0.3 U , and an amplitude a « 0(1.0). For low 
speed air experiments, this criterion indicates that the drag reduction effect 
can only be observed in compliant walls made of elastomers and carefully designed 
to obtain the desired wavemotions under given flow conditions. Indeed, most 
of the reported "successful" experiments were in compliant walls using thin 
membranes . 

According to the proposed model, for a flow speed of 21.3 m/sec (70 ft/sec), 
the structural motion favorable to a compliant effect has a Wavelength of 
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1.68 ram (0.066 in.), an amplitude of 16.76 pm (660 pin.), and a frequency 
of 3700 Hz. 

THE EXPERIMENTAL PROGRAM AT LANGLEY RESEARCH CENTER 

The compliant wall experiments conducted by the Fluid Mechanics Branch, 
High Speed Aircraft Division of Langley Research Center can be classified into 
two categories of test conditions: (1) Wind tunnel test conducted in a Low 

Turbulence Pressure Tunnel (LTPT) with flow speeds of 64 m/sec (210 ft/sec) 
and higher. The results of these experiments have been reported (ref. 5). 

(2) Tests conducted in a 7 in. by 11 in. low speed tunnel where flow speeds 
can be varied between 15 and 45 m/sec (50 to 150 ft/sec). 

Because of the differences in test conditions, model geometries are scaled 
differently for these experiments. So far, LTPT tests have been made in lamin- 
ated structures only and the 7 in. by 11 in. tunnel tests have been done on a 
wide variety of models. 

Currently, drag measurements are made using a direct drag balance, and 
facilities for monitoring surface motions during an actual run are available. 
Thus, in some cases, comparison between the theoretical data presented in this 
report and the experimental values has been possible. 

OVERALL PERSPECTIVE OF THE STRUCTURAL ANALYSIS PROGRAM 
RELATED TO COMPLIANT WALL EXPERIMENTS 


Material Requirements 

The compliant surfaces used for testing have been built of elastomers. 
The compliant wall materials used in the Langley experiments were carefully 
tested in the Materials Division of Langley Research Center. Materials were 
tested in film foam on a tensile testing machine. Typical compliant materials 
used in the experiments were PVC plastisol, Alathon, latex rubber, aluminized 
PVC, polyethylene and aluminized mylar. A brief summary of material testing 
is reported elsewhere (see Appendix A) . 
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Structural Response Analysis 

Previous investigations of compliant wall structural interaction (refs. 7, 

8, 9) have been based on speculation that the compliance of the structure acts 
as a stabilizer of the transitional instability. Neither the effect of supports 
nor the presence of nonlinear effects in the motion received any attention in 
these works. In some later works (ref. 10) the viscoelastic nature of the material 
has been included in the analysis. However, the problem of fluid-structural 
interaction has largely been given only cursory attention. 

The present work has attempted to remedy this situation in two directions: 

1. Modelling of the actual structure has been done carefully; wherever 
numerical simulations are appropriate they have been used. When semi-analytic 
methods are useful we have employed these. Identification of the various loads 
acting on the structure under test conditions is another major problem given 
attention in this report. 

2. Inclusion of structural nonlinearities or other gross effects which 
can affect the actual motion has been made after proper evaluation (for 
instance, if the loading is moderately small, structural nonlinearity need not 
be included when the nonlinear effects are mainly geometric) . If acoustic 
effects are important these effects need to be included in the analysis. If 
mean flow effects come into the picture (flutter, divergence) then the appro- 
priate loading due to this needs to be included. 

Because of the wide variety of compliant models, we resorted to a classi- 
fication which covers basically most of the compliant experiments reported in 
the literature: 

1. Laminated structures: Usually membranes bonded to substrates which 

are rigidly fixed at the bottom. The membranes may be under tension; 

2. Membfanes under deep cavities: The cavities could be fluid filled; 

3. Membranes under narrow air gaps; and 

4. Periodic structures: The periodicity due to multiplicity of equally 

spaced supports usually in the flow direction. 
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STRUCTURAL OSCILLATIONS IN THE PRESENCE OF FLOW 


The response of a structure to loading in the presence of a fluid is quite 
different from its in vacuo response. The motion of the structure introduces 
a disturbance field in the fluid which itself is felt by the structure. Thus 
there are two distinct boundary value problems that need to be solved: 

1. The boundary value problem of the structure of volume V, with the 

given dynamic loads external to it and body forces: only a part of its total 

boundary surface S --say, Sj --interacts with the fluid. 

2. A boundary value problem for a fluid region R interacting with the 
elastic body, the solution of which has the form 

p = p(U, x) Cl) 

where p is the hydrodynamic pressure at the point X on and U is the 

displacement of the points of the body on the boundary. 

We thus break up the total loading acting on the structure into three 
components : 

a. P e , the external pressure loading acting on the structure causing 
the initial or primary motion. 

b. P 1 , the induced pressure loading on the structure due to its inter- 
action with the primary fluid. 

c. pk, the back pressure or the induced pressure from the back side of 
of the structure. 


External Force Field 

For compliant wall experiments under consideration in this report, the 

0 

external pressure loading P can be broken down into three components: 

1. A random convected pressure load due to the turbulent boundary layer. 
(The pressure pulses travel at speeds between 0.5 to 0.8 times free stream 
velocity,) Bull (ref. 11) has made measurements of the spectral pattern of 
this pressure for subsonic flows. At lowest wave numbers, Bull's spectrum 
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does not drop off, and ve have used the dropoff .suggested by. Von Harman and 
Lin (ref. 12), i.e., E(k,t) *» F(I,e°,k,t) where e is the eddy diffusivity, 
and I is Loitsianskii's integral. 

Wherever only qualitative trends are desired (e.g., expected value of 
surface amplitude, power spectral response) we have used Bull's spectrum (with 
the proper modifications) for analysis. 

Recently, we have developed a full simulation of the random pressure field 
using Monte Carlo techniques (ref. 13) . In some analyses using numerical simu- 
lation for structural response (specifically membranes under thin air gaps, 
laminated structures using NASTRAN) , we have used the simulated pressure field. 
In figures 1 and 2 we display the modified wall pressure spectrum and the simu- 
lated wall pressure history at points for a low speed turbulent boundary layer. 

2. Static pressure differentials. In some wind tunnels there will be a 
difference between the static pressure of the fluid and the ambient pressure, 
which will be a function of the flow speed. A knowledge of this deviation is 
essential since it can cause primary surface deformations. 

3. Pressure gradients. The presence of a static pressure gradient in the 
tunnel affects the theoretical predictions in two ways. In all of the analyses 
presented here, a tacit assumption is made that the mean flow is parallel shear 
flow. The nonparallel effect due to gradients requires a nonparallel stability 
model for the proposed model (ref. 1) . Furthermore, the induced pressure force 
evaluation also suffers from lack of consideration of the nonparallelness. 
Hence, if pressure gradients are present in a given experiment, we shall ignore 
that experiment altogether. 

Induced Pressure Field P 1 

The motion of the structure introduces a perturbation pressure field into 
the flow and hence on itself. Sophisticated fluid theories can be developed to 
evaluate the boundary value problem of a flow of infinite extent--one boundary 
of which changes with time. Since many of the experiments were conducted in 
the predivergence regime we have accepted a lower order theory, viz. piston 
theory (potential flow model) , for developing the expression for induced 
pressure on the structure. 


6 






The inclusion of the induced pressure field is unwarranted in some cases. 
These are cases where the divergence speed of the model is far above the opera- 
ting speed. An analysis which neglected these forces would still yield accurate 
results for structural response in these cases. As we proceed in these analy- 
ses we shall indicate the inclusion (or otherwise) for different experiments 
justifying our reasons for same. 

b 

Back Pressure P 

Some of the models tested contained fluid filled cavities beneath them. We 
shall develop theoretical expressions for evaluating these in the next sections. 

ANALYSIS OF A MEMBRANE OVER A DEEP CAVITY WITH EXTERNAL 
FLOW OVER THE MEMBRANE 

We consider (fig. 3) a membrane of flexural rigidity D, tension T, 
simply supported on ends x = 0, x = a, and y = 0 and y = b, and of thickness 
h and density £. 



Figure 3. Membrane over a deep cavity. 
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The dynamic equation of motion for this membrane is given by 

ph -^SL + c f DV 2 V 2 w - TV 2 w = P e + P 1 + P b (2) 

3t 2 3t 


In equation (2) we have neglected effects of midplane stretching. When non- 
linear (stretching) effects need to be taken into consideration, we use the 
following equations: 


ph 


3^w „ 3w 


C * DV^w 


3t 


2 9t - V x 9x 2 


-( N > 


♦ N ^ * N = P e * P 1 * P b 

y 3y 2 xy 3x3y/ 


(3) 


oh — * c i- * MM l— * M-*) ^ 

3t 2 3t Vl - V 2 / I 3x 2 V 2 / 8y2 \ 2 / 


3 2 v 

3x3y 


and 


3w 


3x 


r/3 2 W 

/i - v\ 

3 2 w\ , 1 + v 

3 2 w 

3w ) 

\3x 2 + 

V 2 ) 

M 2 

3x3y 

3 y) 

3v / 

Eh \ 

j 3 2 v , /I - v> 

i 3 2 v 

i 

4* 

3t 


lay* V 2 / 

1 3x 2 

T 


2 3x3y 


(4) 


and 


3w / 3 2 w 1 - v 3 2 w \ 1 + v 3 2 w 3w I 

+ 3y \3y 2 2 3x 2 / 2 3x3yax J 

"-■(*) 1 ( 9 -- 5 ) •#)'••©']! 

/ 

■-■( 99 ) 1 ( 5 -- 9 ) •€©••(©! 


(5) 


(6a) 


(6b) 


N _ Eh 1 3u + 3v_ + 3w 3w 
xy 2(1 + v ) | 3y 3x 3x 3y 


(6c) 
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Linear Model 


We shall assume that the surface motion is made of normal modes in space 
and has the form 


w = £ £ w ft) sin si 

< Z— t rrm J a 


irnrx . mry 
sin — r^- 


Furthermore, we assume that a harmonic dependence on time exists, i.e.. 


W ft) = a e 
mn v mn 


For the linear problem at hand, equation (8) is a justifiable assumption 
as long as P e can be broken down into its harmonic components, i.e.. 


iw.t 

P e =IM 3 


Evaluation of the Back Pressure 

Under the assumption that the fluid motion in the cavity is a perturbed 
state of motion and is irrotational , we define expansions for the velocity 
potential and the elevation of the surface Si as 

<f> = e + e 2 + ... 

, - z<°> ♦ z z™ + ... ) 

<j>^) satisfies the wave equation 


y**Ck) 


c z at" 


where c is the speed of sound in the medium. 
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The boundary conditions are, at n> i.e., on 

[n = z - w(x,y,t) ] 

dn _ ~ _ 8w 8 4> 8w 3 8w 8(f) 

dt 8t 8x 8x 8y 8y " 8z 


and at z = H 


8(f) 

8n 


0 , 


i.e 


9 


34 > 

3z 


0 


( 11 ) 


( 12 ) 


Applying the expansion, equation (9) into equation (11) 


3z 


( 0 ) 


at 


= o 




3t 


3z 


yielding 


*C°> = constant = 0 , 


. . ( 0 ) ( 1 ) 

and since z '+ez =w , 


3w _ 3 4> | 
3t = 3z ' 


(13) 


z=0 


Let us now assume that the solution of equation (10) with the boundary condi- 
tions [eqs. (11) and (12)] has the form 


, . > ^ . nvrrx . nTty 

<l> = a(z) & mn (t) sm — sin -g- 


(14) 
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Thus 


■M - {(=)* (t) 2 | - 4 *W 


with a' (0) = iui 
a’(-H) = 0 


( 15 ) 


Putting 


mn 


■ i@) * (if) 2 - 


0)‘ 


the solution for <J> is given by 


cosh[x (z - H)] 

= T /.I == 


. imrx . n7ry 

= 10) — r rr^- a ft) Sin Sin 

X sinhfX H) mn^ J a b 

mn v mn J 


(16) 


In order to evaluate the back pressure, P, we apply Bernoulli* s equation 
at z = 0, 


+ m = 3_ 

9t p. 


gw 


(17) 


or 


| co coth X H | 

Pb ■ * X a m 1 Sln — Sln b 

I mn I 


(18) 


Evaluation of the Induced Pressure Due to the Primary Flow 
For the velocity potential we can write the governing equation as 
2 


tGht* u » It) ♦ - v2 + 


(19) 


13 


From momentum balance in z direction for the perturbed Eulerian flow 


p ■ -P& * ». £}* 


( 20 ) 


The boundary conditions on 4> are 


|fL = 0 > lfl zss0 = (w + u ~£) w » 


( 21 ) 


using equation (20) we rewrite equations (19) and (21) as 




(22) 


with 


9 p i _ r 9 ti 9 1 2 

9z'z=0 pa L9t 9xJ W * 


(23) 


|E-| - o 

9z 1 z=°° 


(24) 


Solution for P 1 

Once again we use for w the form given by equations (7) and (8) . Because 
of the nature of the boundary conditions [eqs. (23) and (24)] and the governing 
equation (22) itself, we can seek the solution p at the wall z = 0 for a 
mode a) mn (call it P mn ) an d using superposition obtain 

pl - £ E P mn (25) 

m n 


We shall therefore consider the expression for w as 


„ iu)t . mirx . niry 

w = a e sin — — sin - tt - 

mn a b 


(26) 


14 



and rewrite 


mn p i 

w = — [wi + W 2 + W3 + W4J 

and we write 
a 

w = -4— [Wj + w 2 + W 3 + W4] 

where 

i(k)[x cos 0. + y sin 0. + cit] 

= e 3 3 , j = 1,2, 3, 4 

and 0j's are such that 

0 ! = -0 

02 = IT - 0 

0 3 = 0 

04 = IT + 0 


( 27 ) 


( 28 ) 


( 29 ) 


( 30 ) 


where 


0 = Min 



( 31 ) 


k| = 22- g 

a mn 


( 32 ) 


mn 




and 


ci = 


o> toa 


k| mir 3 


mn 


( 33 ) 


( 34 ) 
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In order to evaluate the induced pressure at the wall we compute the pressure 
due to w^ . 

The solution of equation (2) with equations (23) and (24) for an input 
wall function w is given by 


P(z,x,y,t) =£ £ £ 
m n j = l 


mn 

a. 

J 


-|k|p.z i|k|[x cos 9. + y sin 9 . - cjt] 
5 J e J J 


(35) 


where 



(36) 


mn 

a. 

J 


-pjk|[ci + U cos 9j ] 2 


V /Ci + U cos 9 .\ 2 

“(■ F-^) 


a . 


mn 


j 4 


(37) 


with aj =1, a 2 = 1, 03 = -1, and 04 = -1. The wall induced pressure is, 
therefore, given by 


4 


m n j=l 


1 k [x cos 0 . + y sm 0 . + citl 
mn 1 1 1 j j 

a. e 
] 


(38) 


Evaluation of the Divergence Speeds for Subsonic 
Flows Where M << 1 

We shall investigate the nature of the induced pressure P 1 for sub- 
sonic flows with M « 1: 

1. For a) = 0 (for m = 1 , n = 1) we obtain 


i 

P 


a x 


mn 



k 


M 2 

A - M 2 


sin 


1H7TX 

a 


sin 


mr y 
b 


(39) 
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2. For nonzero a) 


p = P , + i Pt 
r real imaginary 

The imaginary term acts as a severe damping term. 

The critical speed can therefore be defined as the speed at which the 
membrane goes into zero frequency flutter. Since frequency is zero, the 
lowest critical speed for the structure occurs when the induced pressure p 
is equal to the stiffness of the structure. Thus, 

° W « (f ) 2 ] 2 - t [(if * gj] m2 

»a‘ 2 [(i) 2 * ( s ) 2 ] V2 

For a given ratio of the flutter boundary (i.e., critical speed) can 

be found from equation (40) . 


Nonlinear Effects 

The linear model is only capable of specifying the divergence speed. In 
an actual experiment one would see membrane divergence and small frequency 
flutter around the critical speed. This is because the nonlinear effects 
have not been included in developing equation (38); also, the critical speed 
will be slightly higher than that predicted by equation (38) since the effect 
of nonlinearity is to increase the stiffness. In order to obtain a valid 
analysis one has to: 

1. Assume that a mn are slowly varying functions of time, and conduct 
an analysis using Kelvin's stationary phase method to obtain closed form 
solutions; or 

2. Develop numerical analytic solutions from p based on Fourier 
transform techniques (see ref. 14). 
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STRUCTURAL ANALYSIS FOR LAMINATED STRUCTURES 

In many compliant experiments reported in the literature (refs. 4 and 5), 
the compliant model consisted of a very thin membrane (thickness ~ 25 ym 
bonded to a soft foam (E ~ 1.4 x io 6 N/m 2 ) of thickness 6 to 12 mm. The 
analysis of the "laminated" structure is performed using proper structural 
modelling. In the following we present different approaches in modelling 
the structure and their basic merits. 


Approximate Analytical Model 

Consider (fig. 4) a membrane (a x b x h) with properties (P^, E m> v ) 
under uniform tension T which is bonded to a substrate (a x b x H) of prop- 
erties (P . E , v ) . We assume that the substrate can be modelled as a semi- 
v s s s 

infinite foundation undergoing plain strain deformations to external loading 
p(x,y,t). The end effects near the edges are thus neglected. The effective 
foundation properties are 



(41) 


Of the load p(x,y,t) acting on the membrane we assume that a portion q(x,y,t) 
is transmitted to the foundation. The deformation state of the foundation is 
therefore assumed as 


u(x,y,z,t) = v(x,y,z,t) = 0.0 
w(x,y,z,t) = W(x,y,t)iJ;(z) 


(42) 


where W(x,y,t) is the deformation of the midplane of the membrane. 
The function (z) is such that 


<|/( 0 ) = 1.0 
ifi (H) =0.0 


(43) 
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The functional form chosen for \jj is 


iKz) 


sinh u(H - y) 
sinh uH 


( 44 ) 


where u is a foundation constant. The stress components of the foundation 
are 



= E f 8w 
T zy 2(1 + v^) 8y 


E f 8W 

2(1 + v f ) * 8y 


(45) 


- E f 8w 
T zx 2(1 + v^) 8x 


2(1 + 


v £ ) 


, 8W 

* 87 


Applying the principle of virtual work, the equation of motion of the substrate 
under the given loading q(x,y,t) is obtained as 


■f T 


a^'dz - 


f 

r\ 


8 Z W 

8t 2 


^ 2 dz = -q(x,y,t) (46) 


Equation (46) can be written as 


q(x,y,t) = m* + k*W - t*v 2 w 
8t 2 


(47) 
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where 


m* = p, 


r H 

. I ty 2 d2 
•Jo 


k* = 


(1 


“ — 2 — f ” t *’] 2 dz > 

- vt) •'o [ 


T* = 


_f H 

V *.) J 


2(1 + v f ) „ Q 


i|) 2 dz 


( 48 ) 


The equation of motion of the membrane can now be developed as 


p h 2-OL + dv 4 W - TV 2 W = p(x,y,t) - q(x,y,t) 
m 9t 2 


(49) 


where 


D = flexural rigidity = 


E h 3 
m 


12(1 - v ) 
v nr 


(50) 


or, using equation (48) , 


(p h + m*) + DV 4 W - (T + T*)V 2 W + k*W = P(x,y,t) 


m 


9t^ 


(51) 


The eigenvalues of the structure can be determined from equation (51) as 


0) 


R 

M * (t) 2 ] 

+ (T + T*)| 

!(") 2 *0?) 2 ] 


mn 


(p m h + m*) 


(52) 
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One of the chief merits of an analysis such as the one above is its simplicity. 
However, the fidelity of the solution so developed depends strongly on the 
correct choice of the foundation parameter u. 

Models Based on Simulation Techniques 

An implicit assumption in the preceding section was the one-dimensional 
(depthwise) variation of all deformations. The foundation thus acts basically 
as a series of vertical springs supporting the membrane. The finite nature 
of the structure makes the approximation carried out in that section rather 
inaccurate. A proper approach to modelling would be using a numerical 
technique, viz. finite elements. 

Using NASTRAN (acronym for NASA Stru ctural Analysis Program) we per- 
formed simulation studies on the laminated structure. The basic models 
are briefly described below. 

Two-dimensional model . The structure is assumed to be two-dimensional 
(X-Z) in the direction of the flow (fig. 5a). The membrane is modelled as 
bar elements (CBAR) with only flexural properties. The substrate is assumed 
to develop only shear deformations due to the external loading; hence, quad- 
rilateral membrane elements (CQDMEM) are used to model that. OFFSET cards 
are used to account for the offset of a node of the membrane element of the 
substrate from the node of the bar element at the interface. 

Plate-spring model . A static three-dimensional analysis is made for an 
elastic slab (substrate) using three-dimensional isoparametric elements 
(CIHEX3) (see fig. 5b) . From the static analysis the equivalent spring 
stiffnesses of the foundation are obtained. The membrane is modelled as 
a plate with membrane action (CQUAD1) . For a given loading the structure 
(CQUAD1-CELAS) is analyzed using NASTRAN. The inertia effects of the sub- 
strate are lumped to grid points using C0NM2 cards. 

Fully three-dimensional model . The substrate is made up of three- 
dimensional isoparametric elements (NASTRAN level 15-9) CHEXA2, and the 
membrane is modelled by CQUAD1 . The offset between the plate grid points 
and the surface grid point of the three-dimensional elements is neglected. 
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Figure 5 




CQUAD1 


30 29 28 



CIHEX3 


(b) 


Quadrilateral plate element and isoparametric element 
with 32 grid points. 



MEMBRANES UNDER NARROW AIR GAPS 


One of the major considerations in developing compliant models for tests 
in the new 7 x 11 wind tunnel has been the ability of the surfaces to undergo 
small wavelength motions. Such motions are "favorable" for drag reduction 
according to the criterion of the Bushnell mechanism (ref. 1) . In order to 
facilitate small wavelength motions, a membrane under narrow a narrow air gap 
was tested at the Langley 7 x 11 tunnel. A drag reduction of about 10 per- 
cent was obtained in this test. The compliant model is shown in figure 6a. 

In figure 6b we indicate the development of small wavelengths due to chopping 
at the bottom. In the absence of the narrow air gap the deflected profile 
would have been as indicated in figure 6c. 

The governing equation of motion for a membrane under a narrow air gap is 
2 

ph — + DV 2 V 2 w - TV 2 w - P(x,y,t) (53) 

8t 2 


with 

w(x,y,t) = -6 

U,y,t) - 0 J t x <_ t < t 2 

where 

(w (x , y , 1 1 ) + s) <_ 0 

and 


(54) 


t 2 - ti + e 

where e is a very short duration of time whose value can be obtained 
through analysis of a dynamic Hertz contact problem. 

The motion of the membrane under these conditions is highly nonlinear. 
Hence no close form solutions are possible. A finite difference solution 
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MEMBRANE (a x b x h, E,p,v) 




DEFLECTED 

SHAPE 


(b) Development of small wavelengths 



(c) Deflected profile in air (for an 
unrestrained membrane) 


Figure 6. Membrane over narrow air gap. 



for equation (53) is possible, and an algorithm with leapfrog differencing 
in time and center differencing in space was used to solve equation (53) . 
The C-F-L limitation imposes 


Ax 

At 


! C sh 


( 55 ) 


For analysis, the bending rigidity of the membrane is neglected. An appro- 
priate choice for At has been made as 


At 



where 


C 


sh 



(56) 


For analysis Ax was chosen equal to Ay (uniform mesh spacing) , and the value of 
Ax was chosen to obtain resolution of the smallest, wavelength possible. The numeri- 
cal model consisted of a 0.25 m x 0.25 m (ID in. x 10 in.) membrane with an air gap 
thickness, 5 = 2.54 ym (10 -lt in.). The smallest wavelength we wanted to capture was 
of length X - 50 mm (2 in.). Thus, the mesh spacing was set at Ax = 5 mm (0.2 in.). 

A one-dimensional analysis of the same problem is possible with greater 
resolution of wavelengths. Here one would assume variations in the y direc- 
tion to be unimportant. However, the one-dimensional analysis can only show 
trends, since the nonlinear nature of the problem prohibits comparisons 
between the one-dimensional analogue of equations (53) and (54) , and the 
two-dimensional problem. 

The conclusions to be drawn from extensive simulation studies conducted 
are that (1) the model equations (53) and (54) have unique solutions (con- 
vergence with reduction of step size) and (2) the solution for dynamic 
response exhibits small wavelength, high frequency motions compared to a 
pure membrane. 
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PERIODICALLY SUPPORTED STRUCTURES 


Compliant drag reduction has been reported in the literature for period- 
ically supported structures (ref. 2) . The length of the bays for these com- 
pliant surfaces corresponds to the wavelengths predicted by the Bushnell 
mechanism. We develop below an analytic method of solution for periodic 
structures based on the work of Mead and Pujara (ref. 15). 

The theory developed here is basically a linear small deflection theory 
with the requirement from consideration that the midplane strains are negli- 
gible, i.e., ^ << 1. We will therefore discard solutions which give ampli- 
tudes w ~ 0(h) or greater. The above requirement is essential, as the ratio 
of the linear restoring force versus the nonlinear stretching force needs 
to be very large, i.e.. 


w x 


Eh 3 


12(1 - v 2 ) 


(¥)' 


>> w 3 Eh 


(¥)‘ 


or 


w 2 

Z— « 
h 2 


1 


(57) 


We also assume a priori that the structure is one-dimensional, i.e., crossflow 
directional variations are neglected (a/b << 1) . 

Figure 7a shows a periodically supported beam, over which is a flow with 
speed U ro . Figure 7b shows an equivalent representation for flexibility at 
the supports and sitting on an elastic foundation (stiffnes = k^) . If the 
foundation is viscoelastic the representation for the foundation is made 
k^ = k^(l + i n f ) where n£ is a typical "loss factor" for the foundation 
and k^ is the complex stiffness of the foundation. 

We shall develop the theory of response of the structure for two cases: 

1. The foundation is absent; instead there is a cavity of fluid at the 
backing, and 

2. The case with the foundation. 
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Case 1. To develop expressions for this case one has to evaluate the 
reaction of the cavity on the bay under consideration. The approach 

is very similar to evaluation of the induced pressure p acting on the 
side where there is flow and hence will not be dwelled upon in detail here. 
Also, one has to set the foundation reaction (k f = 0) .to zero. 


Case 2. Consider an external force field 
x 


-1U 


P = P e 
o 


Z iut 
e 


( 58 ) 


We seek harmonic solutions to the excitation. The transverse displacement 
can be written as 


w(x) 



A 


n 


-i(p+2mr)|- 

e 


icot 


e 


(59) 


We shall assume that the flexural rigidity of the beam is D' 


D’ = D(1 + in b ) 


(60) 


Eh 3 


where n is the loss factor for the beam and D = - 

12(1 - v 2 ) 

of motion of the beam under consideration of excitation is 


The equation 


D* 


— + p, h — + k' w = P e (x,t) + P x (x,t) + Ri 6 (x) + R 2 <$(x - &) 
3x 4 b 9t 2 f 


(61) 


where Ri and R 2 are the reactions of the supports and P 1 is the induced 
pressure. The induced pressure can be evaluated as was done in previous 
sections as 


p 1 (x,t) 



P ito 2 
ik 

n 


A e 
n 


i[cot-un(j)] 


(62) 


where pj is the density of the fluid medium. 
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( 63 ) 



and 'a' equals the speed of sound in the medium, 

y n = y Q + 2mr (64) 


In order to obtain solution of equation (61) for given loading [equation 

i(ym j - ict) 

(58)], we evaluate the virtual work done by virtual displacement 6A m e 
and equate the sum of all virtual work to zero. Because of the periodicity 
of the bays the total virtual work done by the supports to be included in a 
single bay is just due to that of one of the supports. 


The virtual work done by the support at x = 0 is 


6 “r 


(5 Am 


i 




A + K 


E 



(65) 


The virtual work done by the other reactive forces can be obtained by multi- 
plying equation (61) by the virtual displacement and integrating over the 
interval 0 to £. Finally, we indicate the set of simultaneous equations 
obtained by setting the sum of all virtual work to zero as 


m 



P b h + 



or 




A 


n 


+ 



E A „ 


y y 

M m^n 


0 


when m ^ 0 


= p when m = 0 
r o 


( 66 ) 
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The simultaneous set of equation (66) can be solved finally to obtain 
A^ . In practice, only a few of these A^'s have to be evaluated. Taking 
n to be between -6 to +6 is sufficient to obtain the desired accuracy. 
Usually the rotational rigidity of the supports is zero (simply supported, 
i.e., Kr = 0) we indicate here the approach which is simple to solve 
equation (66) for this case. 

We rewrite equation (66) as 

a <t> + ^- y; a = p 

m m l n o 
= 0 

where 



m 


From the case m = 0 we obtain 


m f 0 
m = 0 


(67) 
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P /<#> 
o Y o 
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Then 


-k'K . 

m &<j> * 

T m 



K'k 


for the case where K^, •> » equation (69) yields in the limit 



(69) 


(70) 
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The evaluation of A i leads to the knowledge of w(x) for a given harmonic 
excitation. 

In order to obtain the solution for a wide band excitation one must carry 
out a Fourier analysis of the exciting field and sum up, using the principle 
of superposition, the harmonic responses given by equation (59) . The wide 
band response can then be analyzed to obtain power spectral data of expected 
response. The mean square response is defined as 

E[ojoj*] = e|p(co) P* (co) j H(jco)H*(jw) (71) 

where H(jco) is the frequency response function and E | P (co) P* (co) | is the 
power spectrum of the turbulent pressure field. 

NUMERICAL RESULTS 

In the preceding sections we have developed appropriate techniques for 
analysis of a generic of compliant structures. Using these analyses we 
examine in this section the nature of structural response under reported 
experimental conditions . 


Laminated Structures 

The compliant model tested in the LTPT tunnel at Langley consisted of a 
1.3 x 0.46 m 2 (51 in. x 18 in;), model surface. The compliant surface was made 
of 25 pm (1 mil) mylar [Young's modulus — 3.8 x io -8 kg/m 2 (5 x io 5 psi) , v = 0.3] 
with a backing substrate which is a soft foam (80 PPI polyurethane foam, 

E - 160000 kg/m 2 (200 psi) of thickness 6.25 mm (0.25 in.). In table 1 we report 
the eigenvalues (natural frequency of vibration) obtained using different models 
discussed in the Structural Analysis for Laminated Structures section. These 
eigenvalues are arranged in ascending order. Those eigenvalues listed in the 
last column under "Membrane" represent the natural frequency of the membrane (in 
vacuo) if the substrate were absent. Figure 8 presents the frequency response 
of the laminated structure and indicates that the amplitudes are small; the 
turbulent pressure excites the structure into a wide band response with most 
of the power included in the frequency range below its lowest natural frequency. 
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Table 1. Natural frequencies of a laminated structure.* 


Eigenvalues 

Using NASTRAN 
Models 

Using Eq. (52) 

Membrane 
In Vacuo 

(1,1) 

484.6 Hz 

535.2 

125.6 

(2,1) 

485.0 

536.0 

145.0 

(3,1) 

485.7 

537.6 

172.5 

(4,1) 

486.6 

539.5 

204.9 

(5,1) 

486.6 

541.8 

240.2 

(1,2) 

487.7 

541.8 

240.5 

(2,2) 

488.9 

542.6 

251.2 

(3,2) 

490.4 

543.9 

268.0 

(4,2) 

490.5 

544.6 

277.4 

(6,1) 

492.0 

544.9 

290.0 


* The membrane is 25.4 pm (1 mil) mylar, substrate 6.25 mm (0.25 in.) thick, 
80 PPI PU foam £E = 1600 kg/m 2 (2.2 psi) , v = 0.1]. The model is 1.3 m 
x 0.46 m ( 51 in. x 18 in.) and with a tension on the membrane of 40 kg/m 
(2 lbs/in.). 

** The foundation parameter u is chosen to be 2.0 
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We conclude that laminated structures do not show any promise as "passive" 
compliant surfaces on the basis of the Bushnell mechanism. 

Membranes Under Narrow Air Gap 

As was indicated in the Membranes Under Narrow Air Gaps section above, 
the analysis of the membrane over a narrow air gap only presents qualitative 
information regarding the nature of the motion. Using the simulation program 
for pressure (ref. 13) we have developed solutions for a 0.25 m x 0.25 m 
compliant model with 25 ym mylar over 2.5 ym air gap and a flow of 
15 m/sec. Such flow conditions typically occur in the 7 in. x 11 in. tunnel 
at Langley. The boundary layer thickness was approximately 12 mm. The 
pressure history and displacement for the midpoint of the surface is 
indicated in figure 9 for a duration of time. 

Figure 10 indicates the influence of spatial discretization on the solu- 
tion obtained. Figure 11 presents an analysis of the frequency response of 
the motion. It is shown that the narrow air gap-membrane configuration 
yields high frequency motion. 

Membranes Under Large Cavities 

In some reported experiments (ref. 3) membranes over fluid-filled cavities 
have shown drag reduction. In order to understand the nature of the surface 
motion we examined a recent compliant wall experiment: the membrane is 25 ym 

thick mylar and the cavity is 6 mm deep and filled with air. The test 
speed is 15 m/sec. The mean square response with frequency is indicated 
in figure 12. It is shown that the motion is purely low frequency, large 
wavelength. Therefore, we do not believe that membranes over large air gaps 
are promising candidates for compliant experiments. 

Figure 13 indicates the response of a membrane with a water-filled 
cavity. Again the presence of the cavity shifts the response curve towards 
the low frequency end of the spectrum. 
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Figure 9. The pressure history and resulting surface motion 
for a membrane over an air gap (U = 15.2 m/sec, 
25 mm boundary layer thickness) . 
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Figure 10. Influence of resolution on predicted surface motion. 
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Figure 12. Membrane over large air gap — frequency response curve 







Periodic Structures 


The original reporting of compliant drag reduction (ref. 2) was for a 
periodic structure made of a rubber diaphragm (2 mm thick) supported on 
periodic stubs of rubber with water filling between the stubs. The height 
of the stubs was 1 nun and the width 0.5 mm. Using the above data we 
examined the surface response to turbulent flow (medium is water) at a speed 
of 15m/sec corresponding to the experiment of Kramer. Figure 14 indicates 
the surface response. We determined that the periodic Kramer surface was a 
high frequency passive surface and could have been a drag-reducing surface 
satisfying the conditions of Bushnell criterion. In nondimens ional units 
our analysis indicates for the Kramer surface a A + = 800, h + = 100, and 
c = 0.7 U . The analysis includes a nominal value of damping for the 
structure. We expect, however, that the actual damping in the structure 
may limit the amplitude to h + — 20 or 30 in the experiment. Such an exci- 
tation is capable of reducing the rate of burst production and hence could 
have been responsible for the observed drag reduction. 

Most attempts toward reproducing the results of Kramer over the years 
have been unsuccessful . One of the chief reasons for this has been the 
belief that Kramer's experiment delays transition; hence, subsequent experi- 
ments at various laboratories were conducted as transition experiments. In 
a few cases the flow had been accidentally tripped to turbulence, and there 
have been reports that drag reduction was indeed observed for these cases. 
However, the majority of data points showed no favorable drag changes (ref. 
16) . We also emphasize that proper attention was not given to the structural 
motions required. Consequently, the geometries tested were not facsimiles 
or scale models of the Kramer surface. There is, therefore, a need for re- 
examining Kramer's surface for possible drag reduction. 

PASSIVE WALLS FOR 7 IN. x H IN. TUNNEL 

The 7 in. x 11 in. tunnel at Langley has a speed range of 15 to 45 
m/sec. The flow is tripped about 0.3 m from the test model using rougheners 
(sandpaper) . Boundary layer surveys on the test model indicate that a rela- 
tively thick boundary layer (6 = 6 to 18 mm) is formed on the model. 
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Figure 14. Steady-state response of Kramer stubbed peri 
(U = 18 m/sec water flow) . 
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The model itself is 0,275 x 0.9 m and. floats on a drag balance (see fig. 15). 

The facing plates smoothly fan the flow onto the test area. When rigid plate 
measurements are made the whole surface area is available. However, the 
compliant section on which measurements are made spans only a 0.175 x 0.4 m 
area in order to avoid nonuniformities of flow around the side walls. The 
adjustable backplate provides a facility for adjusting the cavity depth which 
can be filled by a foam or with any fluid. When a backing such as air is 
used, the back pressure is controlled using a vacuum pump. The overall 
pressure gradient in the test region is kept at a minimum. However, there 
is a slight pressure gradient of the order of 1 to 5 N/m 3 in the tunnel test 
section. The total pressure on the flow side fluctuates slightly from the ambient, 
and the tunnel static pressure differs from the ambient pressure as a small 
function of the flow speed under conditions of optimal operation. While 
evaluating the response of structures under flow we included this small 
static pressure differential acting on the structure in addition to the 
existing force fields. 

The vacuum section chamber in the model (fig. 15) is used to put tension 
on the compliant wall. The widths of the cavity of the suction chamber are 
adjusted in such a way that the tension is uniform. This gives the ratio 
of the cavity widths in x and y directions as 



(see Appendix B) in order to obtain uniform tension of the membrane. By 
directly observing the center deflection of the membrane in the cavity on 
any side and noticing the pressure differential, one can directly evaluate 
the tension in the membrane knowing the properties of the membrane. 

Finally, direct calibration can be made between the tension in the mem- 
brane and the vacuum pressure provided that the properties of the compliant 
membrane are fully known. 

In figure 16 we show the observed surface motion using a schlieren system 
for taking area photographs. The surface is a membrane over a narrow air 
gap. The dark spots in the picture correspond to points where the membrane 
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Figure 16. Area photograph of the response of a membrane 
over a narrow air gap; backing surface is 
10 PPI foam. 
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touches the 10 PPI foam surface after the excursion through the narrow air 
gap. The arrows indicate the direction of flow. The tension in the membrane 
is low. 

In figure 17 we show the surface motion as a function of time at a point 
on the surface as seen in the oscilloscope. Each division in the screen 
corresponds to 5 milliseconds (in the x direction) and in the ordinate 
1 division corresponds to 0.1 volts of photodetector output for the bottom 
picture and 0.01 volts for the top one. The output of the photodetector is 
directly proportional to the maximum angles of the surface motion and has 
been calibrated already. 

From figure 17 we learn that the surface motion of the membrane over the 
narrow air gap is high frequency (between 400 Hz and 2 KHz) and has ampli- 
tudes of 2 to 3 yms. 

In figure 18 we show the area photograph of a ’’dulcimer" model. The 
model consists of a membrane sitting on periodically supporting strings. There 
is a separation between the membrane and the strings * and during surface 
oscillation the membrane slaps on the strings. The dark lines correspond 
to the strings in contact with the membrane. The strings are kept in tension. 
The motion of this surface is again small wavelength and at moderate fre- 
quencies. Figure 19 shows the motion of a point with respect to time. The 
time scale is 20 milliseconds/division. The flow speed is U^ = 48 m/sec. 

Again, the surface amplitudes are only about 5 yms (1 volt»2.5 ym) . 

Using the theoretical analyses in the preceding section, we also con- 
ducted some experiments where the wall was made to flutter. (The spacings 
between supports were designed to be such that low frequency flutter occurred.) 
One of the chief problems in these experiments was adjusting the pressure in 
the back chamber continually while the experiment was going on, in order to 
suppress the panel divergence. 

The models were tested for drag reduction too. In most of these compli- 
ant experiments there was little drag change. The observed motions did not 
suggest that the surface motion was in the range of parameters suggested by 
the Bushnell mechanism either. Further compliant experiments where compli- 
ant motions suggested by the Bushnell mechanism can be produced (active walls/ 
passive walls) are underway. 
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Figure 17. Surface motion of a point on the membrane over 
narrow air gap as function of time. 
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Figure 18. Area photograph of the surface motion of a 
"dulcimer” model. 
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Figure 19. Surface response with time for the dulcimer model. 
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CONCLUSIONS 


Theoretical analyses for common compliant wall geometries are presented. 
Analyses of surface motion for the compliant experiments indicate that the 
original Kramer experiments conducted in water could have interrupted the 
turbulent burst mechanism and produced the drag reduction. Membranes over 
thin air gaps can also be used to produce the small wavelength, large fre- 
quency motion required for air experiments. It is possible that in some 
laminated surface compliant wall experiments, delamination of the structure 
occurred and thereby created a narrow air gap situation; thus the observed 
drag reduction in the experiment of Walters (ref. 4) and the Langley experi- 
ment (ref. 5). The membranes over large air gaps do not show any promise 
at all. The "observed" drag reduction reported in reference 3 can be due to 
other effects, as suggested in reference 17. It seems rather difficult to 
obtain using passive walls the desired surface motion for most low speed and 
transonic speed experiments in air. Thus it may well be worthwhile to use 
active walls, i.e., walls where the motion is driven onto the surface for 
compliant experiments. Finally, there is a clear need to understand the 
entire fabric of fluid structural interaction, and this can only be achieved 
by bringing in and coalescing structural expertise with the fluid mechanics 
of turbulence. We believe that this report goes as far as making a first 
tentative step in the right direction toward such a goal. 
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APPENDIX A 


FILM MATERIAL TESTS 


By 

Todd Hodges 1 


INTRODUCTION 

Tests have been completed on eight materials in film form to support the 
compliant skin drag reduction project in the Fluid Mechanics Branch. The 
Fluid Mechanics Branch was interested in the elastic modulus of the materials 
they were using in their wind tunnel models . 

TEST EQUIPMENT AND MATERIALS 

Tests were conducted on the tensile testing machine located in room 003 of 
building 1293A. The C range load cell was used. Depending on the material 
being tested, the full scale ranges used were 4.5 N, 9N and 22 N (1 lb, 2 lb, 
and 5 lb). Accuracy of the load cell is equal to, or exceeds, 0.25 percent of 
the range in use or 0.50 percent of the indicated load, whichever is greater. 
The amplifier and recorder accuracy is within 0.50 percent of full scale range. 

Air grips were used to hold the test specimens. Test gage length was 
75 mm (3 in.), and specimen width was 12.5 mm (1/2 in.). Film thickness 
was measured by a power-driven film micrometer. 

The test specimens were mounted and pulled at a constant crosshead speed 
of 12.5 mm/min. The tensile modulus (E) was calculated according to the 
attached equation derivation with an electronic calculator. Five samples 
were used on each test run except for the latex rubber, where only three 
samples survived the film cutter. 


1 Research Associate, Virginia Polytechnic Institute and State University, 
Blacksburg, VA. 
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DISCUSSION 


The following modulus values were measured on the eight systems: 

1. Ethylene vinyl acetate copolymer (Alathon) (tested in the longitud- 
inal direction). E = 138 x 10 8 N/m 2 with a high of 1.42 x 10 8 N/m 2 

J average 6 

and a low of 13.7 x 10 8 N/m 2 . This material showed good consistency in 
thickness and modulus. 

2. Alathon (tested in the transverse direction). E = 107 x 10 8 

v average 

N/m 2 with a high of 114 x 10 8 N/m 2 and a low of 9.04 x 10 7 N/m 2 , also showed 
good consistency in thickness but not as good in modulus. 

3. Latex rubber, E = 2.04 x 10 6 N/m 2 , high 2.16 x 10 6 N/m 2 low 

average 6 

1.88 x io 6 N/m 2 . Lowest modulus of all the materials. We were able to salvage 
three test samples from the film-cutting device. 

4. Aluminized kapton, E = 2.87 x 10 9 N/m 2 , high 3.23 x 10 9 N/m 2 

average 

low = 2.69 x io 9 N/m 2 . Modulus consistency was not good. 

5. Aluminized P.V.C,E =2.5x 10 9 N/m 2 high 2.87 x 10 9 N/m 2 low 

2.04 x io 9 N/m 2 . Thickness and modulus variation were high. 

6. P.V.C., E = 2.3 x io 8 N/m 2 high 2.55 x 10 8 N/m 2 low 1.95 x io 8 

average 6 

N/m 2 . Thickness and modulus variation were high. 

7. Polyethylene, E avera g e = 15 x 10 8 N/m 2 high 1.57 x io 8 N/m 2 low = 1.44 x 
10 8 N/m 2 . Thickness and modulus consistency were good. 

8. Aluminized mylar, E avera g e = 5.11 x 10 9 N/m 2 high 3.7 x 10 9 N/m 2 low 
2.53 x io 9 N/m 2 . Thickness and modulus consistency were not good. 

The aluminized films generally showed less consistency in thickness and 
modulus than plain films. Thic could be due to the coating technique they 
use and its effect on the substrate film. 

Ranking the materials in order of modulus we have: 

1. Latex rubber: E = 2.04 x 10 6 N/m 2 (296 psi) 

2. Alathon (transverse): E &v = 107 x io 8 N/m 2 (15500 psi) 

3. Alathon (longitudinal) : E av = 138 x io 8 N/m 2 (20000 psi) 

4. Polyethylene: E &v = 15 x 10 8 N/m 2 (21700 psi) 
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5. P.V.C. : E = 2.3 

av 

6. Aluminized P.V.C. : 

7. Aluminized leap ton: 

8. Aluminized mylar: 


x 10 8 N/m 2 (33400 psi) 

E av = 2.5 x io 9 N/m 2 (3.6 x 

E = 2.87 x 10 9 N/m 2 (4.16 
av v 

E Q = 3.11 X 10 9 N/m 2 (4.51 
av 


10 5 psi) 
x 10 5 psi) 
10 5 psi) 
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APPENDIX B 


VACUUM SUCTION TECHNIQUE FOR TENSIONING COMPLIANT MEMBRANES 


In this section we shall briefly discuss the vacuum suction technique 
for tensioning 'the compliant membranes. The cross section of the vacuum 
suction chamber is shown in figure Bl. The shaded area of (X x * L ) is the 


D 


C 



SUCTION 

CHAMBER 


Figure Bl . 


working area of the model on which a uniform tension needs to be applied. In 
order to apply uniform tension in the working model the width of the cavities 
j i x , Zy have to be related in some fashion to the dimensions L^, Ly. 

Initially the membrane is taped to the sides ABCD with care such that 
there is little tension on the membrane. The suction chamber pressure is at 
ambient pressure during this time. Now the chamber pressure is adjusted to 
be lower than the ambient pressure by a magnitude Ap. The deflections of 
the membrane at the center of the channels wj, W 2 are noted. 


In figure B2 we indicate schematically what happens to a cross section 
of the membrane. For the simply supported membrane of length z and width 


fL + 2Z ) where the ratio 
v y y J 

equation as: 


/ \ 

\K + 2 V 


<< 1 we can approximate the governing 
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Figure B2. 


T — = Ap 
3x 2 


(B-l) 


w(0) = w(Jl x ) = 0 


(B-2) 


Assuming 


" = E w m ^in t 


mirx 


(B-3) 


one obtains for .W m , after substituting in equation (B-l) 


W = 
m 


2Ap 
(mir) ' 


(l - C-D m ) 


(B-4) 


Thus all odd modes are present in the solution and 


Wi = T W sin 
1 m 


imrx 


m=l 


4Ap£^ 

tt 3 T 


L 

m=l,3,5 . . .» 


imrx \ 1 

sin T "I ~T 

Z x / m 3 


(0.125) \=rj Ap 


(B-5) 
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Similarly 


w 2 


0.125 



Ap 


Thus 


Wi 

«2 




(B-6) 


(B-7) 


Now, the stretching that the membrane undergoes in the cavity is con- 

tributed by the stretching of the membrane in the working area. The stretch- 
ing in the cavity SL , is given by 



2 2 


(B-8) 
x 


= 0 . 0417Ap 


/A 3 x 

fe 


16Ap£ 


ii4t x 2 


E - 


nr 


The stretching AJ, is assumed to be caused by a uniform loading P per unit 
width applied to the membrane of the working area of the model, in which case 

PL 2 

-=£- = A£ = 0.0417Ap 
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or 


P = 




(B-9) 


Similarly the loading on the end in the y direction is 

p ‘=(f) °- 0417 (^) 


CB-10) 


Finally, we would expect to have P = P* and T = T for uniform tensioning 

x y 

when 



L vl/3 


© 


(B-ll) 


which is as 


given in equation (72 ) . 
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APPENDIX C 


CODE FOR ANALYZING MEMBRANE WITH BACKING, 
OR WITH CAVITIES BEHIND IT 




PROGRAM MA I N ( I NPUT * OUTPUT . T APE5 = 1 NPUT • TAPE6=0UTPUT ) 

COMPLEX DU.X2.X3 

D I MENS ION AK ( 75 <24 > i AF ACT ( 75 ) .BF ACT ( 24 ) 

DIMENSION W 1(5500) 

DIMENSION OMC ( 75 . 24 ) . CM1 (75.24) 

DIMENSION FF ( 80 ) 

C*********** ************************* ******************************************* 

C 

c 

c 

C THIS PROGRAM EVALUATES THE RESPONSE OF A MEMBRANE WITH DEEP 

1 CAVITY! AIR) « A MEMBRANE WITH WATER BACKING. AND AMEMBRANE 

2 ON ELASTIC FOUNDATION 

1 CA V I TY (AIR). AMEBRANE 

C BULL-S SPECTRUM IS USED TO OBTAIN THE EXCITATION FIELD 

C** INDUCED AERODYNAMIC EFFECTS USING PISTON THEORY. 

C** EFFECTS OF BACKING USING POTENTIAL FLOW THEORY. 

C** THEORETICAL MODEL USES LINEAR THEORY 

C** PROPER VALUES OF DAMPING INCLUDED USING LAB ESTIMATES 

C OUTPUT FREQUENCY. MEAN SQUARE FORCING FUNCTION .MEAN SQ 

1UARE RESPONSE FUNCTION. ME ANSQARE AMPLITUDE. 

C ALL DATA TO BE IN IN-LB-SEC UNITS 


i 



c 

c 

c 

NR = 5SNW=6SG = 386»0SEPS = 0« 00000000 1 SC 0=1 . 35* 1 0 . **4SP I =3 . 1 4 1 5926 
DO 110 J=1 .5 

READ (NR .111 1 ( FF ( I >.1=1*80) 

1 1 1 F ORMAT ( 80 A 1 ) 

110 CONTINUE 

WRITE (NW* 088 ) 

808 FORMA! < //OX.* CASE 1 WTI ' IKK Ml M. MEM Of* 0 A V . C.3 JML'M PLASH) 

DO 777 I JKL= l .3 

READ (NR. 101 ) A * B . T * EM 1 . ETA . DELT .UREF . U ♦ RHO 1 , RH02.H 

101 FORMAT (8F1 0.5 ) 

READ ( NR . 1 02 ) I UK . AKF 

102 FORMATt 1 1 .F14.6) 

EMI =EM1 * ( . 1 )**7*RH01 =RH0l / (12. ) **4£RH02=RH02/ (12. )**4 

0STAR=7./72.*DELT 

DO 5 L= 1 . 75 

AL1 = (L— 1 ) *2 • + 1 . 

DO 5 N= 1 .24 

AN 1 = ( N— 1 . ) *2 • + 1 » 


AK1 = <AL1*P1/A)**2+(AN1*PI/B > **2 



0MC(L.N)=C0*SQRT(AK1 ) 


AK CL.N )=AKI 

5 CONTINUE 

DO 6 L= 1 .75.2 
OMCUT=OMC ( 1 * 1 ) — 5 « 

6 AFACT(L)=I»0 
DO 7 L = 2 . 74 » 2 

7 AFACT ( L ) = — 1 • 0 
DO 8 N= 1 .24.2 

8 BFAC.T ( N ) = 1 • 0 
DO 9 N = 2 « 24 . 2 

9 BFACT < N ) = — 1 .0 
WRITE(NW.80l )U 

801 FORMAT (//20X.F15. 6 ) 

DST 1 =DSTAR 
Q=0.5*RH01 *U**2 
C0NST=Q**2*DST1/U*1 . 

CONST1 =0.0064516 
CONI =DST1/U 
UC=0.8*U 

802 FORMAT ( //20X . *l.OWER CUTOFF FREO. . (BULL ) * . F 1 5 . 6 .*UPPFR CUTOFF * .F 1 5. 



D0M=1 .0$0M2=2.*Pl *DOM 
0MEG=0.05/C0N1 

EE1 =EXP <0. 1 )$EE2=EXP (0.0235 ) $EE3=EXP (0.40 ) 
DO 2 J=1 ,5 
0M3=0M2 
0M2=0M3*1 0. 

DO 31=1 .90 
JK= JK+ 1 
A t = 1 

OMl =0M3* (1 . +A 1/10. ) 

0ML=0. 1 *0M3 
F=OMl/2./PI 

OMC1 =250.*1 . /CON 1 *2 • *P I 
OMCC=OMC 1/10. 

IF (OMl .GT.OMCC )GO TO 28 
IF (OMl .GT.OMC1 ) GO TO 31 
IF (OMl .LE.OMEG )GO TO 23 
AT0M=8. *0M1 *C0N1 
TW0M = 2 . *OM 1 * CON 1 
FRACOM=0.47*OM1 *C0NI1 
IF (A TOM— 12. ) 12. 13,13 


El =1 ./EXP (ATOM ) 



GO TO 14 


13 E1=0.0 

14 CONTINUE 

I F ( TWOM- 1 2 • ) 15.16.16 

15 E2=l ./EXP(TWOM) 

GO TO 17 

16 E2 = 0 • 0 

17 CONTINUE 

I F ( FRA COM- 1 2 . 1 18, 19. 19 

18 E3=l ./EXP (FRACOM) 

GO TO 20 

19 E3=0 . P 

20 CONTINUE 
E=E2*3.7+0.8*E3~3.4*E 1 
GO TO 22 

28 E= 3. 15/225.0 
GO TO 22 
31 E=0 • 0 

GO TO 22 

23 E= <0Ml/0MEG)**4*(3.7/EEl-3.4/EE3+0.8/EE2 ) 

22 CONTINUE 

FORE=CONST*E 

803 FORMAT (//20X.*F0RC I NG FIELD — AMP . = * , F 1 5 . 6 , *FREQ. =* , F 1 5. 6 ) 



BK = 0M1/UC$BK 1 = 1 ,+COS (BK*A ) $BK2=S I N ( BK+A ) 
DU=CMPLX 10.0,0.0) 

ETAl = ET A *0M 1 
OMN = 10M1/C0 )**2 

XX = PI /A1XXX=XX**2$ASQ=A**2S6KSQ=BK**2 
IFUJK.EO.l J PH02 = 0.0 
DO 10 L=1 ,75 
AL= (L-l • )*2.+l . 

ALl =AL*XX 
APU=ABS (Al 1 -BK ) 

IF (APL.LE.EPS 1G0 TO 51 
ALX1 = t AL 1 **2-BKS0 )*ASO 
D 1 1 = RK 1 *8. *AL. * AF ACT (L ) / AI.X 1 
D 1 2 = D1 1 *BK2/BK] 

GO TO 52 
D1 1 =0. 

D12=-4./PI*AFACT(l_) 

CONTINUE 
DO 10 N=1 ,24 
AN=(N-1 • ) * 2 • + 1 , 

XI 1 =D 1 1/AN*BFACT(N) 

X12=D12/AN*BFACT IN) 


X3=CMPLX(X1 1 ,X12 ) 



Al< 1 = AK ( L » N ) -OMN 


IF (OM1 .LT.OMCUT )GO TO 9.3 
IF(ABS(AK1 >.LE.EPS)GO TO 92 
I F ( AK 1 -EPS ) 9 1 i 92. 93 

91 AK2=-AK1 

AK2 = SQRT ( A:<2 >*H 

CT = COS ( AK2 )/S I N < AK2 ) 

XI =-<EMl+RH02*H*CT/AKP )*OM1**2+AK(L.N)*T+RH02*G 
XI =X1 +AKF 

X5*ETA1+0MI **2*RH0 1 *H* 1 ./AK2 
X2 = CMPI_X (XI »X5 ) 

GO TO 97 

93 AK2 = SQRT(AK1 )*H*2» 

CX=£XP< AK2 ) 

CTH= ( CX+ 1 . )/<CX-l • ) 

XI =-( (RHOl +RH02*CTH ) *2 • *H/AK2+EM 1 )*OM 1 **2+A< <L . N ) * T +G*RH02 

XI =X1 +AKF 

X2*CMPLX (XI » ETA 1 ) 

97 CONTINUE 

CM 1 (l_*N)=CABS ( X3/X2 ) 

DU=X3/X2+DU 

o. 

92 CONTINUE 
10 CONTINUE 



987 FORMAT ( 1BE8.2) 

Wl ( JK >=CABS(DU)**2 

^ W1 ( JK+1 000 ) = Wi ( JK )*FORE 
00 

Wl (JK+1 000 ) = CONST 1 *W1 ( JK+1 00 0 ) 

W= W+W 1 (JK+1 000 )*OML 

WRITE (NW* 909 >F .FORE . Wl ( JK > . Wl (JK+1 000 > . W 
909 FORMAT (BX.5E20.7 > 

3 CONTINUE 
2 CONTINUE 
1 CONTINUE 
777 CONTINUE 
STOP 
END 

C GIVE EMI =EM1*1 0,**7$RH0S =RH0*144.**2 

DATA TO BE FED ARE LENGTH OF MODEL < A ) . W 1 DTH OF MODEL ( B )♦ TENS I ON OF MEMBRANE (T ) 
MASS OF MEMBRANEfEMl ) .DAMPING! VI SCOUS ) COEFF 1 C I ENT ( ETA ). BOUND .L AVER THKNSS(DELT) 
.REF.VEL (UREF ) . VELOC ITY (U) .DENSI TY FLO 1 ( RH01 1 . DENS .FLD2 ( RH02 ) HT OF CAVITY(H) 
FORMAT NUMBER IS 101 .8F10.4 


26.0 

8.0 

0.0?. 

1.1 

0.0002 

0,6 

636.0 

636.0 

0,00234 

1 #9379 

0.5 






u\ 

o 

• 

o 

18.0 

1 #0 

1 • 1 

0.0002 

1 .6063 

1318.9 

1318.9 



0.00234 0. COP 34 0.00 

50.0 18.0 1.0 1.1 

0.00234 0.00234 0.08 

1 0.666666 


as 

so 


0.0002 1.6063 1310.0 1318.9 




APPENDIX D 


CODE FOR ANALYZING PERIODIC STRUCTURE RESPONSE 




PROGRAM KRAMER ( INPUT, OUTPUT , TAPE5= I NPUT , T APE6 = 0UTPUT , TAPE 1 ) 


C* ************ **R****************************** ***************** **************** 
C** 

C THIS PROGRAM EVALUATES THE MEAN SQUARE RESPONSE OF A KRAMER 

1 PERIODIC SURFACE. THE SURFACE IS A RUBBER DIAPHRAGM SUPPORTED 

2 PERIODICALLY (80 MIL APART) BY RUBBER STUBS AND FILLED WITHBWATER. 

C** INDUCED AERODYNAMIC EFFECTS USING PISTON THEORY. 

C** EFFECTS OF BACKING USING POTENTIAL FLOW THEORY. 

C** THEORETICAL MODEL USES LINEAR THEORY 

C** PROPER VALUES OF DAMPING INCLUDED USING LAB ESTIMATES 

C ANALYSIS --SEE SECTION 6. 

C ALL UNITS LB-SEC-FT 

C BULLS SPECTRUM FOR EXCITING FIELD 

C** 

C***************R**************** *********************** ******* ************ ***** 
DIMENSION A (21 ) *B <21 ) «C <21 ) ,D(21 ) . 1 J<21 > » FORE < 1 2000 ) , W ( 1 5 ) * 

1 A SUM (15) 

DIMENSION GE < 2 1 ) , F ( 2 1 ) 

DIMENSION HEF (21 ) 

NR=5$NW=6 
C DATA 

w H=0,040*RHO=62»4$G=3B6.*U=720.4ETA=0.00l iE=2000.SEMl =8 . 8S ANU =0. 45 



DELTA= 1 .08SY = 0.3>TP I =6. 283185260*00 = 57000. 

ETA=0. 0000001 
EPS=0« 0000000 1 

C INITIALISATION of data 

EM1=EM1*< .1 )**6S’RH0 = RH0/32.2/l2./12./12./12.$DSTl=7./72.*DELTA 
AKF = E* A 1 =H/AKF $Pt = A 1 
AKF=20.0SA1 =H/AKF&Pl =A1 
EMI =EM 1 *A 1 

DD=E*H**3/3./ < 1 • — ANU**2 ) 

D1 =DD*A1*ETA1 =ETA*A 1 *UC = 0 . 8*U*0=0 . 5*RH0*U**2*C0N 1 =DSTl/U$01 =Q**2*C 
1 ONI 

G 1 =G*RH0*A1 

0MEG=D. 05/C0N 1 SD0M=0. 1 *TP I S0M2=D0M 

DOM= 1 0 • *TP I 

0M2=D0M 

OMCl =25 » *TP I *U/DELT A 
OMCC=10.*OMC1 
C INITIALISE ARRAYS 

DO 101 J= 1 ♦ 15 
ASUM ( J ) =OSW ( J ) =0. 

101 CONTINUE 

C MAJOR LOOP FOR INCRENENT OF OM 


JJ=0$EE1 = EXP (0. 1 0 )SEE2=EXP ( 0 .0235 ) *EE3=EXP ( 0 . 4 0 ) 



C DO 2 JK = 1 ,6 

0M3=0M2$0M2=0M3*1 0. 

C DO 1 JL=1 .90 

DO 1 JL=1, 20000 
A J = JL 

0M1=0M3*< 1 .+AJ/10. ) 

FF = 0M 1 /TP I 
JJ= JJ+1 
UC=0.8*U 

1 F < OM 1 • LE • OMEG ) GO TO 23 
IF (OM1 .GE.OMC1 )GO TO 511 
IF (OM1 .GT.OMCC >GO TO 513 
AT0M=8.0*0M1*C0N1 
TW0M=2 • *OMl *C0N1 
FRAC0M=0.47*0M1 *CONl 
1 F ( ATOM— 60 • 1 12, 13,13 

12 E 1 = 1 • /EXP (ATOM ) 

GO TO 14 

13 El =0. 

14 CONTINUE 

IF (TW0M-60* >15,16,16 

15 E2=l ./EXP (TWOM ) 


GO TO 17 



16 E2 = 0.0 

17 CONTINUE 

IF (FRACOM-60. ) 18,19. 19 

18 E3=l ./EXP (FRACOM) 

GO TO 20 

19 E3=0. 

20 CONTINUE 
E=F2*3.7+0.8*E3-3.4*E1 
GO TO 22 

511 E=3. 18/225. 

UC=0 .39*U 
GO TO 22 
513 STOP 

23 E=3.7/EE1 +0.8/EE2-3.4/EE3 

E=E* (OM1 /OMEG )**4 
22 CONTINUE 

FORE (UJ ) =01 *E*0. 0000 1 

AKt = (0M1 /CO )**2S>AMU = 0MI/UC*H$ETA I =ETA 1 #0M 1 
M=-ll SSUM=O.SSUMl =0. 

DO 30 Ml =1 .21 
M=M1 -1 1 
AMI =M 


AML= ( AMU + AM1 *TPI )/H 



AK=AML#*2-AKl 


A ( Ml )=D1*AML**4+G1-0M1**2*EM1$B(M1 )=ETAI 
IF(ABS(AK).LE.EPS)GO TO 40 
IF(AK-EPS>39»40*41 

40 I J (M 1 )=1 
GO TO 47 

39 CONTINUE 
I J (Ml )=0 
AK=-AK 
AK=SQRT (AK) 

S«=AK*H 

Si =SIN(SK )$C1 =COS(SK ) 

IF (SI .EQ.O ) GO TO 901 

A ( Ml ) =A (Ml )+RH0*A1*C1/S1/AK*0M1**2 $B(Mi )=B(M1 )+Rh 0*A 1 /AK*OMl **2 
GO TO 902 

901 I J ( M 1 >=1 

902 CONTINUE 
GO TO 47 

41 I J (Ml >=0 
AK=SQRT(AK) 

SK=AK*2.*H 
SK 1 =EXP (SK ) 

CTH= (SKI +1 . )/ (SKI -1 . ) 


I 



oo 


52 CONTINUE 

A(M1 )=A(Ml )-0Ml**2*RH0/AK*Al*(l ,+CTH) 

47 CONTINUE 

IF ( I J ( M 1 ) ) 48 *48 * 49 

48 CONTINUE 

A2=A (Ml )**2+B (Ml )**2 
SUM=SUM+A (Ml )/A2 
SUM 1 =SUM1 +B (Ml )/A2 

49 CONTINUE 
30 CONTINUE 

A9 = A ( 1 1 )**2+B (11) **2 

GAM 1 = A 1 * < A ( I 1 )* (SUM+1 . ) +0 ( 1 1 )*SnM] ) 

GAM 1 =GAM1 /A9/ ( ( 1 ,+SUM >**2 + SUM 1**2) 

GAM2 = A 1 / A9* ( A ( 1 1 )*SUM1 -B ( 1 1 >* (SUM+i , ) )/ (SUM1 **2+ ( 1 .+SUM>**2 > 

DO 60 Ml =1 ,21 

IF ( I J(M1 ) )61 *62,61 

61 CONTINUE 
D ( M 1 ) = 0.0 
C ( M I ) = 0 . 0 
GO TO 67 

62 CONTINUE 

C (Ml )= ( A (Ml )*GAMi +B (Ml )*GAM2 )/ i AIM] )**•*; + tM M 1 , -> , 

D ( M 1 ) = ( B ( M 1 ) *GAM 1 -A ( Ml ) *GAM2 ) / ( A < Ml ) **2+B ( M 1 > **2 ) 



92 CONTINUE 

WR ITE(NW»200 )FF • ASUM (4 ) » GE U 1 > » GE < l 2 ) *GE < 1 0 > • HEF ( 1 1 ) » HEF ( 12 ) iHEFU 
10 ) ,W(4 ) ,F0RE( JJ ) 

603 FORMAT(/50X»*IM* ) 

209 FORMAT ( 1 X « 1 0E1 2 *5 ) 

208 FORMAT (//5X. 10E12. 6 ) 

1 CONTINUE 

2 CONTINUE j 

201 FORMAT <1X< 7E 16*7) I 

. j 

STOP 

END 


^4 

KD 




APPENDIX E 


CODE FOR ANALYZING MEMBRANE UNDER NARROW AIR GAP 




■ -’•' 33 ! 


JOB, 1 ,2500.77000430000. A4677 R4623 100718 BIN34 

RUN(S ) 

REQUEST, TAPE1 ,HY. SAVTP , R I L . RS6 , TU/20 • 

REWIND! TAPE 1 1 

PROGRAM MEMGAP ( 1 NPUT , OUTPUT , TAPE5= 1 NPUT. TAPE6 = 0UTPUT , TAPE 1 ) 
DIMENSION W<51 .51 ) ,W0<51 ,'il ) , WN ( 5 1 . 5 1 ) , P ( 5 1 , 5 1 ) ,AU< 101 ) 

C T = TENS I ON (LBF/IN) 

C RHO= MEMBRANE DENSITY (LBM/CU.ET) 

C E = SUBSTRATE MODULUS OF ELASTICITY c PS I ) 

C D=SunSTRATE DEPTH (IN) 

C GAP=GAP BF TWEEN MEMBRANE AND SUBSTRATE (IN) 

C H=MEMBRANE THICKNESS (IN) 

C B=DAMP I NG COEFFICIENT ( LBF-SEC/CU . FT ) 

C DX=ELEMENT SIZE (IN) 

C NI=NUMBER OF NODES IN X-OIREC.TION 

C NJ = NUMBER OF NODES I N. Y-D I RECT J ON 

C NT = Nl.JMBER OF TIME STEPS — OT=DX/C WHERE C IS WAVE SPEED 

C GM IS TIME STEP CONTROL DT=DX/(C*GM) 

T = 1 

RH0=87 

00 

w 

F=i . 


D=n.?5 



GAP=0 • 000 1 


NW = 6 

H = 0 » 00 1 

B=0. 000001 

NI = B1 SNJsBl $DX=0.2 

(jsBOI-DFLT A = 1 • f’P 1 = 3 • 1 41B92G 

NT=20000 

PMAX=0.00234*U**2*0.5*0.0l 
PI = Pf/AX*U**2/DX**2 
PI =P1 * ] 44 . 

F = 1 .?*U*12./PI/2./DELTA 
c convert input TO CONSISTF-nT units 

TS 12. *T 

RHO=RHO/32. 174 
E=144.*£ 

D=D/12 

GAP=GAP/12. 

H=H/ 1 2 . 

C WAVE SPEED IS GIVEN BY SO. RT. OF T/(RHO*H) 

C= (T / CRH0*H ) )**0.B 
DX=DX/1 2 
nr=F')fo . i #-*6 


GM=0X/C/DT 



INITIAL I Z" 01 SPl. ACF VI: NTS 


DO '> 1 = 1 iNl 

DO 1 J= 1 iNJ 
WO (I tJ)=0, 

W ( I * J ) = 0 . 

'AN ( I . J ) =0 . 

P ( I . J ) = 0 . 

C.ONT INUE 
CONT I NUE 

CALCHLATF NEW DISPLACEMENTS 
PRF'.'UPE LOADING CA\ Lit. ENTERED HERE 
I T = 2 

N I X = m I - I 
NJV=NJ- 1 

A1 =R*GM/ (2.*T*DX )+ (GM/DX )**2 
A2= 1 ./DX*1 ,/DX 

A3=H*GM/ (2. »T*OX ) - (GM/DX )**2 
A <4 = F / ( D * T ) 

AS=1 ./T 

TE= 1 • /F10M=?.*PI #F$T 1 =DX/U'ETWOT 1 =2.*T1 
TPLT1 =TE+TI 


TPLT2=TF+TW0T1 



r 


EN I =N I — 1 

T M A X = E N I /K *DX+1 ./F' + DX/U 
8 REWIND 1 

DO 130 N=1 .NT 
I T=2+N 
T j mf=N*P r 

T = T I me 

DO 22 7 JX = ?,N'I 
DO 227 JY-2.NJ 
227 P ( I X, JY ) =0. 0 

IF (T IMF-TVAX )20 1 ,210,2 10 
2.01 DO 209 I V - 2 , N I X 
EXX= I X-2 
T = T I ME-EXX* nx/l I 
I F ( T 1 1 52 , 1 , 1 F4 

154 I F ( T-T 1 ) 1 55 , IFF, 151 

155 CONTINUE 

DO 221 JY = 2 , NJ Y 

221 P ( IX, JY >=P1 * ( T/OM-S IN (OM*T )/0M**2 ) 
GO TO 152 

151 IF (T-TWOT 1 ) 1 56, 156, 1 57 

156 CONTINUE 


DO 222 JY=2 , NJY 



222 P(IX.JY) = P1*< TWOT 1 /OW-T/OM + 2 . *S I N (OM* I T-T 1 ) ) /0M**2-. C > I \i ( OM k T ) /0M-)l*2 
1 ) 

GO TO 1?? 

157 IF(T-TE) 158, 158, 159 

158 CONTI NUF 

DO 223 JY=2,NJY 

223 P ( 1 X, JY ) =P1 *2 ,/OM**2*S I N (DM* ( T-T 1 ) ) * ( 1 . -COS < OM*T 1 ) ) 

GO TO 15? 

159 IF (T-TPLT1 ) 160, 160, 1 61 

160 CONTINUE 

DO 224 JY = 2 , NJY 

224 PI IX, JY ) =P1 * ( I T-T 1 ) /OM+2 ,/0M**2 *S I N ( OM* ( T-T 1 ) ) -S IN ( OM* I T-TWOTl ) I/O 
1 M**2 ) 

GO TO 15? 

161 IF (T-TPLT? ) 1 62, 162 , 1 5? 

162 CONTI NUF 

DO 225 JY = 2 , NJY 

225 P( IX, JY )=P1 * < ( T-TE-TWOTl l/OM-SIN (OM* (T-TWOTl ) )/0M**2 ) 

152 CONTINUE 

209 CONTINUE 

210 CONTINUE 

00 

^ DO 20 I =? ,NI X 


DO 10 J = 2 , NJY 



SUMS'.*.' ( l + l 4 J J+W ( !-l • J )+W( N J+l )+W ( I • J-l ) 


SUM=SUM+(2#*GM#*2-4» )*W< 1 i J) 

SUM] =W0( 1 , J > 

SUM2=A4*(W( I ♦ J)+GAP) 

WN< I ♦ J)=SUM*A2/A1+SUM1*AT/A1-P( 1 * J)*A5/Ai 
IF (W ( I • J )+GAP)5«5 *6 
WN( I • J)=2«*W0< I • J)-W( I • J> 

CONTINUE 
CONT I NUE 
C.ONT I NUE 
A I =N I — 1 
BI =1 00. 

DEL = B I /A I 
I DEL =OEL +0.2 
DO | 9 1=1,101 

El =1-1 

J= ( 1-1 )/IDEL+l 

EJ=J-1 

J1 = J + 1 

XIsF I/DEL-EJ 

AU ( I ) =w N < J « 1 1 ) * ( 1 .-X T ) +WN ( J 1 ,11 ) *X I 
OOMT I NiJ c 


WRITE ( 1 )AU 





C SHIFT LOCATIONS OF 

00 91 I =2 .N I X 
f>'.’ M J=f iNJY 
'VC) ( 1 • J ) = W (I,J) 

W ( I i J ) " to N ( I i J ) 

91 CONTINUE 

130 CONTINUE 

STOP 
FNP 


ARRAYS 


oo 

<0 
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